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Abstract — We investigate the Fourier transform of the current 
through a memristor when the applied-voltage frequency is 
smaller than the characteristic memristor frequency, and the 
memristor shows hysteresis in the current-voltage plane. We find 
that when the hysteresis curve is "smooth", the current Fourier 
transform has weights at odd and even harmonics that decay 
rapidly and monotonically with the order of the harmonic; when 
the hysteresis curve is "sharp", the Fourier transform of the 
current is significantly broader, with non-monotonic weights at 
high harmonics. We present a simple model which shows that 
this qualitative change in the Fourier spectrum is solely driven 
by the saturation of memristance during a voltage cycle, and not 
independently by various system parameters such as applied or 
memristor frequencies, and the non-linear dopant drift. 



I. Introduction 

SINCE the first experimental demonstration in 2008 by 
the Williams group at Hewlett-Packard (TJ, over the past 
four years, memristors |2| and memristive systems [3| have 
been extensively investigated on theoretical and experimental 
fronts H). This ongoing research is spread across diverse 
areas such as modeling memristive systems [5|-[9| (includ- 
ing memcapacitors and meminductors) pO) , pi] , developing 
memristor emulators fl2) , (T3J, the potential use of memristors 
for high-density non-volatile memory p4|-|[T6|, and the devel- 
opment of artifical neural networks |17|-|20|. The existence 
of memristive systems (or a memristor) was postulated nearly 
four decades ago based on symmetry arguments [2|, and many 
examples of their realization, including the charge transport 
across ion channels J2T| , were discussed [3|. However, the 
first simple, physical model that leads to a charge- (or flux-) 
dependent resistance was developed by the Williams group (Tl; 
although this model is primarily applicable to a semiconductor 
thin-film device, it has become a starting point for most of 
the recent theoretical studies. The two unusal and remarkable 
properties of a memristor are its hysteretic behavior in the 
current-voltage characteristics i(v) and the qualitative change 
in its hysteresis loop structure that takes place when the 
frequency of the voltage source is much smaller than the 
characteristic frequency of the memristor in the presence of 
the applied voltage Q 

Historically, hysteretic behavior is associated with magnetic 
materials where it has been extensively explored both the- 
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oretically and experimentally. Since the dependence of the 
magnetic induction B on the field intensity H is too complex 
for an analytical treatment, based on approximate models |22) , 
Fourier analysis has been used to study the dependence of 
magnetization M(H) and B(H) on a sinusoidally varying 
field intensity [23|-[25]. A corresponding Fourier analysis of 
the response of basic electrical circuit elements to a sinu- 
soidally varying voltage is trivial; it is nonzero only at the 
frequency of the applied voltage. Similarly, Fourier analysis 
of non-linear, passive elements is straightforward. Memristors, 
with their hysteretic properties and a simple, analytical model 
of their behavior, are ideally suited for Fourier analysis, 
although such properties are barely explored J26). 

In this paper, we present numerical and analytical results for 
the Fourier response of the current i(t) through a single mem- 
ristor with an applied sinusoidal voltage v(t) — vosm(oj a t). 
Our two salient results are as follows: (i) A smooth current 
hysteresis leads to the presence of both even and odd har- 
monics of the applied voltage frequency ui a , and the weight 
of these harmonics decreases monotonically with their degree, 
(ii) When the current hysteresis changes to a sharp, switching 
behavior, the Fourier transform of the current shows a broad, 
non-monotonic structure; in particular, it has maximum weight 
at the second harmonic over a large range of parameters. 

The plan of the paper is as follows. In Sec. [II] we describe 
the memristor model that takes into account the nonlinear 
dopant drift. We present typical numerical results for the 
Fourier analysis of the current and its correlations with the 



system parameters in Sec. Ill In Sec. IV we present a simple 



analytical model of the hysteretic characteristics and show 
that this model qualitatively explains the numerical results. 
We conclude the paper with a brief discussion in Sec. [V] 



II. Description of the memristor model 

We use the model of a memristor as two resistors in 
series [1|: its effective resistance is given by M(w) = 
(w/D)R on + (1 - w/D)R oS , where w(t)/D = x(t) is the 
fractional width of the doped region with < x < 1, D is 
total width of the thin-film device, R on is the resistance of the 
memristor if the entire device is doped, and i? D ff ^ R on is the 
resistance of the undoped device. The motion of the boundary 
between the doped and the undoped regions is approximated 
by the drift velocity of (oxygen vacancy) dopants, and is char- 
acterized by the dopant mobility no- The equations describing 
the time-evolution of a circuit with a single memristor and a 
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voltage source are given by 



i(t) 
dw 
~dl 



G(w(t))v(t) 
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i(t)F 



D 
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where the polarity r) = ±1 denotes that the doped region 
increases (decreases) with an applied positive voltage, and 
the dimensionless window function F(x) is used to model 
the non-linear drift of dopants near the thin-film boundaries, 
x(t) ~ or x(t) ~ 1 |6J-||9(. The memductance G(w(t)) is 



G(w(t)) = 



Con G, 



on^off 



M(w(t)) G on - x{t)AG' 



(3) 



With Gon — -Ron' '-'off — -R ff ' anC ^ ~~ G on — G a S ^> G ff. 

In the following numerical calculations, we use i? n as the unit 
of resistance, i$ — vo/R on as the unit of current, and note 
that the characteristic frequency of the memristor is given by 
uj m — 2tt/jL£iVq/ D 2 . To simplify the numerical analysis we use 
the frequency of the applied voltage u a as the unit of frequency 
(and its period T a = 2n/uj a as the unit of time), instead of the 
conventional practice where the memristor frequency sets the 
scale Jl], |6j, (9). The hysteretic effects for different values 
of = u a /uj m < 1 are investigated by varying vo and thus 
uj m . Typical memristor parameters |Tj are i? Q n ~ 10 3 Ohm, 
D ~ 10 nm, fin ~ 10~ 14 cm 2 V _1 s _1 , and thus imply that 
the memristor frequency is Lu m ~ 50 KHz and the current is 
io = 1 mA for vq ~ 1 V. 

We use the following one-parameter family of window 
functions 



F p (x) = 1 - (2x - l) 2p 



(4) 



which ensures that the dopant drift is symmetrically sup- 
pressed as x — > 0, 1 |6|. Note that Eq.Q can be generalized 
to non-integer values of the parameter p > 1 by using the 
absolute value of the argument, \2x — 1|. Since the window 
function vanishes at the boundaries, F (0) = = F(l), it has 
two fixed points in the time-evolution of the fractional width; 
however, this is not of physical concern for the following 
reason. Suppose the fractional width changes from an initial 
value x\ — 1 — Si to a new value x 2 = \ — after the passage 
of charge Q where S 2 <C <5i -C 1. This process describes the 
approach to the fixed point x = 1. The charge can be estimated 
from Eq.(|2ji by approximating the window function near the 
boundary, F(x = 1 — 8) ps ApS, 



Apfj, D R on \S 2 



(5) 



Eq.{5]) shows that the charge Q required to reach the fixed 
point, S 2 = 0, diverges. Thus, starting from an initial width 
x 6 (0, 1) it is impossible to come across the "terminal state 
problem" |7), (8) in a finite amount of time. An identical 
argument applies to a memristor where the doped region 
decreases with time, thus approaching the x = fixed point. 
Therefore, in the following, we use the phrase "saturation of 
the fractional width" to denote fractional width x(t) arbitrarily 
close to one (zero) without being exactly equal to one (zero). 



III. Fourier analysis: numerical results 

We calculate the (discrete) Fourier transform (FT) of the 
memristor current i(t) that, in turn, is obtained by numerically 
solving Eqs.([T])-(j2]) over many periods T a = 2n/uj a of the 
applied sinusoidal voltage to ensure stability. The discrete FT 
of the current is given by 



1 N 

-Y 



(6) 



where N is the number of points that approximate the con- 
tinuous function i(t) over a single period, tj — (j — I) At = 
(j - l)T a /N, and v k = kuj a for k = -N/2, N/2 - 1 are 
the N possible discrete frequencies; we verify that our results 
do not depend upon N and thus are truly representative of the 
continuum limit, N -» 0, At -» 0. 
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Fig. 1, The left-hand column shows the amplitude of current FT, as a 

function of frequency Vk/ui a , and the right-hand column shows corresponding 
time-domain i(v) curves. As the ratio of applied voltage frequency to the 
memristor frequency Q = ui a /u) m decreases, the hysteresis changes from 
smooth (top row) to sharp (center and bottom row), and the current FT 
changes from monotonic decay (blue circles) to a broad, bumpy structure 
(red squares, black diamonds) with a maximum weight at the second harmonic 
(black diamonds). 

Figure [T] shows typical results of such an analysis. The 
left-hand column shows the amplitude of the FT of current, 
|Z(i/fc)| = fk)\ an d the right-hand column shows corre- 
sponding time-domain i(v) curves plotted in the units of io 
and vq. Note that since the time -integral of the current over 
a single period is zero, the spectral weight at k — vanishes 
in all cases. These results are obtained with R ff/R n = 100, 
x(t = 0) = 0.5, 7] = +1, and window function parameter 
p = 5. When £1 = u) a /u> m > 1, the memristor behaves 
essentially as an ohmic resistor with constant value and the 
current FT is nonzero only at Vk = ±w - 

As decreases (top row, fl = 1/6) a smooth, pinched 
hysteresis loop develops in the i(v) plane. The maximum 
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current is ~ Ron/Rinit ~ 0.02 where R^t = 0.5R on + 
(1 — 0.5)i? o ff is the initial resistance of the memristor. The 
corresponding current FT (filled blue circles) shows finite, 
monotonically and rapidly decreasing weights at har- 

monics Vf. = kuj a for 1 < k < 10. These even and odd higher 
harmonics are genetically present in the regime where the 
memristor shows hysteresis. When fl = 1/7 (center row), the 
i(v) hysteresis loop now has a butterfly shape, and the max- 
imum current is now determined by the significantly smaller 
resistance value R ~ R on <C i?i n it that the memristor reaches 
when x(t) approaches one. The corresponding current FT 
(filled red squares) shows a broad structure with appreciable, 
non-mono tonic weights at high harmonics v k with k < 25. The 
bottom row shows that for Q = 1/8, the i(v) curve exhibits 
a clear switching behavior from a state with high resistance 
R- m it to a state with very low resistance R ~ R on accompanied 
by the saturation of the fractional width x(t). The Fourier 
transform of the current (filled black diamonds) shows that 
a markedly higher weight at the second harmonic develops 
compared to the weight at the first harmonic. Note that since 
the maximum current through the memristor changes by an 
order of magnitude or more as f2 decreases, the "small spectral 
weight" at high harmonics for il = 1/8 (bottom row) is, in 
absolute terms, comparable to the maximum spectral weight 
of the first harmonic for 51 = 1/6 (top row); if the latter is 
experimentally measurable, so will be the former. 

Thus, the Fourier analysis of the memristor current shows 
two main features. The first is that a smooth hysteresis 
curve is accompanied by a current spectrum with monoton- 
ically decreasing weights at both even and odd harmonics. 
The second is that this spectrum changes to a broad, non- 
monotonic, bumpy structure with nonzero weights at very high 
harmonics as the hysteresis curve changes to a butterfly shape 
characterized by two resistance values. The question is: is it 
the frequency £1 = ui a /uj m that drives this change? 

To investigate the connection, if any, between the broad 
spectrum of the current FT and the saturation of the fractional 
width, we obtain |I(z^-)l f° r different window functions while 
keeping all other parameters constant. The right-hand column 
in Fig. [5] shows the time-evolution of x(t) from its initial 
value x(t = 0) = 0.6 for three values of p = {1,2.5,5}; 
the left-hand column shows corresponding amplitude of the 
current FT. When p = 1 the highly non-linear window function 
reduces the dopant drift velocity and therefore the fractional 
width reaches a maximum x(t = T a /2) ~ 0.9 that is well 
below saturation (top row); the corresponding shows 
monotonically decreasing weights at harmonics v% = ku) a for 
2 < k < 10 (filled blue circles). For p — 2.5 the dopant 
drift is higher, the fractional width barely reaches saturation 
at t = T a /2 (center row), and the current FT now shows a 
broad, non-mono tonic, bumpy structure (filled red squares). 
The bottom row shows that when p = 5, the dopant drift 
is essentially constant, the dimensionless width saturates to 
one for a significant fraction of the period T a , and the current 
spectrum now shows maximum weight at the second harmonic 
(filled black diamonds). We emphasize that these results are 
for a single value of the memristor frequency, Q = 1/5. 

These numerical results strongly suggest that saturation of 
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Fig. 2. The left-hand column shows the amplitude of current FT, 11(^)1 as 
a function of frequency u^/ui a and the right-hand column shows the time- 
evolution of the fractional width x(t) over a single period < t/T a < 1 for 
window functions with p = 1 (top row), p = 2.5 (center row), and p = 5 
(bottom row). These typical results suggest that the qualitative change in the 
current FT is triggered by the saturation of the fractional width. 



the fractional width is instrumental to the qualitative change 
in the current spectrum and the emergence of higher-weight 
second harmonic. They show, in particular, that neither the 
frequency il nor the dopant drift dynamics, characterized by 
the window function index p, solely drive the emergence of 
higher harmonics with increasing weights. 

IV. Fourier analysis: theoretical model 

Armed with these insights, we now present a theoretical 
model that qualitatively explains the properties of current 
Fourier transform in terms of time evolution of the memduc- 
tance. In general, the Fourier transform of current i(t) is given 
by the convolution, 



N/2-1 

£ 

j=-N/2 



(7) 



where Q and V denotes Fourier transforms of the memductance 
G(t) and the applied voltage respectively. For a sinusoidal 
voltage, Eq.([7]) becomes 



(8) 



where we have used v k qp w a = ^fc^i- Thus, an analytical 
model for the memductance FT Q(y k ) can shed light on the 
evolution of the current Fourier transform. 

We start with the case when the fractional width x(t) does 
not saturate (top-right panel in Fig. |2|. Then the current can 
be expressed as sum of two, low-order, odd polynomials in 
voltage 



i(t) - P(v(t))6(v'(t)) + Q{v(t))9{-v'(t)), 



(9) 



4 



where 0(x) is the Heaviside function, and the sign of v'(t) = 
dv/dt denotes whether the voltage is increasing with time or 
decreasing. For example, the hysteresis for Q = 1/6 in Fig. |T| 
can be modeled by using P(x) = a\x + a-sx 3 and Q(x) = 
b\x — 63X 3 where aj,6j > and a\ — G; n it < 61, whereas 
when Q ~ 1, P(x) = aix = Q(x). Although the polynomials 
are odd functions of time, due to the presence of the Heaviside 
function, the Fourier transform of Eq.Q contains both even 
and odd harmonics. We note that although such approximation 
works well for the most part, for any degree polynomial, it 
cannot capture the diverging slope of the i(v) curve that occurs 
at the voltage extrema, v = ±vq. For a power-law current 
with a single exponent, i — A n v n , and sinusoidal voltage, the 
Fourier transform is given by 



piecewise constant function, 



T(v k ) = A, 



(20 



(10) 



The binomial coefficient, and therefore the spectral weight 
decreases monotonically with the order of the har- 
monic 1 < k < n; note that I(vjS) = for k > n. It follows 
from Stirling's approximation that the current FT will decay 
rapidly as exp(— k 2 /n) for large k. This exponential nature 
of the decay does not depend upon the power-law exponent 
n, and is applicable to a general polynomial i = P(v). 
Thus, the two-polynomial model qualitatively explains the 
monotonically decreasing weight at even and odd harmonics 
for the current FT. 
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Fig. 3. Time evolution of the fractional width x(t) (top panel) and the 
memductance G(t)/G n (bottom panel) over a single period < t/T a < 1. 
The system parameters are R ff / R on = 100, p = 5, r] = 1, x(t = 0) = 0.7, 
which implies that Ri nl t = 30.7iJ O n- The results show that although the 
fractional width slowly saturates to one over time (top panel), G(t) increases 
sharply from Gj n it = 0.033G on to G on at time t/T a = a (bottom panel). 

When the fractional width x(t) saturates at low frequencies, 
the memductance G(t) increases from its initial value G- m n = 
R^ n l t -C Gon to a maximum value of G on . Figure [3] shows the 
typical time evolution of the fractional width x(t) (top panel) 
and the corresponding memductance G(t) (bottom panel) as 
a function of 17. In general G(t) can be approximated by a 



G init < t/T a < a, 
G{t) = { G on a < t/T a < 1 - a, 
G in it 1 - a < t/T a < 1, 



(11) 



where 0<a<l/2is determined by the system parameters; 
for example, a rj 0.37 when fi = 1/3 and monotonically 
decreases with ft. Its Fourier transform is given by 



Q(v k ) = G on Sk,o - 2a(G on - Gi n i t ) 



sm(27ra£;) 
2irak 



(12) 



The presence of the sinc-function in Q(vk) and the fact that 
G on ^> Gi n it imply that memductance spectral weight 
is non-monotonic and its higher harmonics decay slowly as 
1/fc for large fc. Then it follows from Eq.([8]l that the broad, 
bumpy structure of the current FT arises from the same source. 

V. CONCLUSIONS 

In this paper, we have presented numerical results for the 
Fourier response of the current through a single memristor in 
the presence of a sinusoidal voltage. With simple models, we 
have shown that the qualitative change in the current FT can 
be traced to the qualitative change in the hysteresis loop in 
the i(v) plane. The results presented here are for a positive 
polarity, ?; = +1, where the fractional width saturates to one. 

When r\ = — 1, the memductance of the device decreases 
from Ginit to G g < Ginit- Thus, to obtain a ratio of the 
initial and saturation memductance values comparable to the 
positive polarity case, we will need an initial fractional width 
x(t) ~ 1. This, in turn, implies that the frequency fl required 
for width saturation down to zero is much smaller than the 
corresponding values in the positive polarity case. Therefore, 
non-monotonic current spectrum and higher-weight second 
harmonic response [26 1 in negative polarity memristors require 
significantly different conditions. 

Our approach can be used to analytically investigate other 
phenomena, such as power in a load resistor Rl in series with 



a memristor, P(t) = i (t)Ri, |26|. Since Rl is constant, the 



form the current-voltage hysteresis and the subsequent results 
remain the same. Since V{vk) — ^jli^k — v j)T(vj) where 
^(vk) is the power FT, our results can be used to understand 
its behavior |26|. 
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